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Abstract 

As exemplified by the Kuramoto model, large systems of coupled oscillators may undergo a 
OO ' transition to phase coherence with increasing coupling strength. It is shown that below the critical 

o ; 

^^ ■ coupling strength for this transition such systems may be expected to exhibit 'echo' phenomena: 
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a stimulation by two successive pulses separated by a time interval r leads to the spontaneous 
formation of response pulses at a time r, 2r, 3r..., after the second stimulus pulse. Analysis 
of this phenomenon, as well as illustrative numerical experiments, are presented. The theoretical 



Q . significance and potential uses of echoes in such systems are discussed. 

pj . PACS numbers: 05.45.Xt, 05.45.-a, 89.75.-k 
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Large systems consisting of many coupled oscillators for which the individual 
natural oscillator frequencies are different naturally occur in a wide variety of 
interesting applications. As shown by Kuramoto, such systems can undergo a 
type of dynamical phase transition such that as the coupling strength is raised 
past a critical value, global synchronous collective behavior results. In this paper 
we show that another interesting, potentially useful, behavior of these systems 
also occurs below the critical coupling strength. Namely, we demonstrate that 
these systems exhibit echo phenomena: If a stimulus pulse is applied at time t = 0, 
followed by a second stimulus pulse at time t = t, then pulse echo responses can 
appear at t = 2r, 3r, .... This phenomenon depends on both nonlinearity and 
memory inherent in the oscillator system, the latter being a consequence of the 
continuous spectrum of the linearized system. 

I. INTRODUCTION 

Due to their occurrence in a wide variety of circumstances, systems consisting of a large 
number of coupled oscillators with different natural oscillation frequencies have been the 
subject of much scientific interest [l|, |2[ . Examples where the study of such systems is thought 



liesjal 



to be relevant are synchronous flashing of fireflies [SJ] and chirping of crickets [J], synchronous 
card,ac pacemaker eellsfl. brain function fl. coordinafon of oscUlatory neurons governing 
circadian rhythms in mammals [7|], entrainment of coupled oscillatory chemically reacting 
cells (si, Josephson 
model of Kuramoto 
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unction circuit arrays[9], etc. The globally-coupled, phase-oscillator 



111 ] exemplifies the key generic feature of large systems of coupled 
oscillators. In particular, Kuramoto considered the case where the distribution function 
of oscillator frequencies was monotonically decreasing away from its peak value, and he 
showed that, as the coupling strength K between the oscillators is increased through a 
critical coupling strength Kc, there is a transition to sustained global cooperative behavior. 
In this state {K > Kc) a suitable average over the oscillator population (this average is 
often called the 'order parameter') exhibits steady macroscopic oscillatory behavior. For 



K < Kc a. stimulus may transiently induce macroscopic oscillations, but the amplitude of 
these coherent oscillations (i.e^the magnitude of the order parameter) decays exponentially 
to zero with increasing time[ll|. In the present paper we consider the Kuramoto model in 
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FIG. 1: Illustration of the echo phenomenon. Stimuli at times t = and t = t lead to direct 
system responses which rapidly decay away followed by echo responses that can arise at times 2r, 
3t, .... The 'response' plotted on the vertical axis is the magnitude of the complex valued order 
parameter, Eq. (8). See Sec. IV for details of this computation. 



details). In order for this phenomenon to occur, the system must have two fundamental 
attributes, nonlinearity and memory. Nonlinearity is necessary because the response seen 
in Fig. 1 is not the same as the sum of the responses to each of the individual stimulus 
pulses in the absence of the other pulse (which is simply the decay that occurs immediately 
after the individual stimuli, without the echo). Memory is necessary in the sense that the 
system state after the decay of the second pulse must somehow encode knowledge of the 
previous history even though the global average of the system state, as represented by the 
order parameter, is approximately the same as before the two pulses were applied. 



Echo phenomena of this type, occurring in systems of many oscillators having a spread 
in their natural oscillation frequencies, have been known for a long time. The first example 
was the 'spin echo' discovered in 1950 by Hahn|12l]. where the distribution of frequencies 
resulted from the position dependence of the precession frequency of nuclear magnetic dipoles 
in an inhomogeneous magnetic field. [The spin echo forms the basis for modern magnetic 
resonance imaging (MRI).] Subsequently, echoes for cyclotron orbits of charged particles in 
a magnetic field have been studied for the cases in which the distribution in frequency was 






due to magnetic field inhomogeneity (l3l| . relativistic dependence of the particle mass on its 






energy [iJ], and Doppler shifts of the cyclotron frequency la] . Another notable case is that 
of plasma waves, where the frequency distribution results from the Dop pler shift of the wave 
frequency felt by charged particles with different streaming velocities 16|]. Although echo 
phenomena are well-known in the above settings, they have so far not received attention in 
the context of the Kuramoto model and its many related situations. It is our purpose in the 
present paper to investigate that problem. Two possible motivations for our study of echoes 
in the Kuramoto model are that they provide increased basic understanding of the model 
and also that they may be of potential use as a basis for future diagnostic measurements of 
related systems (see Sec. V). 

In what follows. Sec. II will give a formulation of the model problem that will be analyzed 
in Sec. Ill and numerically simulated in Sec. IV, while Sec. V will provide a discussion of 
the implications of the results obtained. 

II. FORMULATION 

We consider the basic Kuramoto model supplemented by the addition of a 5-correlated 
noise term n{t) and two impulsive stimuli, one at time t = 0, and the other at time t = t, 

N 

dOi/dt = uJi + K/N ^ sin(ej - Oi) - h{ei)A{t) + n{t) , (1) 

i=i 

A(t) = do6it) + diSit - r) , (2) 

{n{t)n{t')) = 2i5{t - t') , (3) 

/^(0)=5^/^„e*"^ K = h*_^ , ho = , (4) 



where /il„ denotes the complex conjugate of /i_„. In the above 6i{t) represents the angular 
phase of oscillator i, where i = 1,2, . . . ,N ^ 1; and uji is the natural frequency of oscillator 
i where we take Ui for different oscillators (i.e., different i) to be distributed according 
to some given, time-independent distribution function g{uj), where g{uj) has an average 
frequency u = J ujg{uj)duj, is symmetric about uj = Co, and monotonically decreases as 
|c<j — u;| increases. 

To motivate the impulsive stimuli term, consider the example of a population of many 
fireflies, and imagine that the stimuli at t = and at t = r are external flashes of light at 
those times, where the constants do and di in Eq. (2) represent the intensity of these flashes. 
We hypothesize that a firefly will be induced by a stimulus flash to move its flashing phase 
toward synchronism with the stimulus flash. Thus a firefly that has just recently flashed will 
reset its phase by retarding it, while a firefly that was close to flashing will advance its phase. 
The amount of advance or retardation is determined by the 'reset function', h{6). Since the 
reset function h{6) depends on properties of the fireflies, we do not specify it further. Let 
91' and 9^ represent the phases of oscillator i just after and just before a stimulus flash at 
t = or t = r. Then we have from Eq. (1) that 

Letting F{9) = f d9/h{9), we obtain 

9i = F'\d, + F{9r)). (6) 

In our subsequent analysis in Sec. Ill, we will for convenience assume that dp is small, in 
which case {9f — 9^) is small, and we can use the approximation, 

9t = 9r + dX9r)- p = 0, 1 . (7) 

Following Kuramoto we introduce the complex valued order parameter R{t), 

1 ^ 

in terms of which Eq. (1) can be rewritten as 

d9i/dt = Ui + {K/N)Im[e-'^^R{t)] - h{9i)A{t) + n{t) . (9) 



In our analysis in Sec. Ill we will take the limit N ^ oo useful for approximating the 
situation where A^ ^ 1. In that limit it is appropriate to describe the system state by a 
continuous distribution function f{6,u,t), where 

j f{e,u;,t)- = l, (10) 

and the fraction of oscillators with angles and natural frequencies in the ranges {6, 9 + d9) 
and {(jj^uo + doo) is f{6,uj,t)g{uj)dujd6/2TT. The conservation of the number of oscillators then 
gives the time evolution equation for f{9,uj,t), 

^ + ^ {/ [^ + KIm{R{t)e-^') - M^)A(t)] } = ^^ , (H) 

W{t) = j dujf\{u,t)g{u), (12) 

where R* denotes the complex conjugate of i?, f{uj,6,t) = 1 for t < 0, and, in writing 
Eq. (12), /i represents the e*^ component of the Fourier expansion of f{uj, 6, t) in 9, 
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f{0O,e,t)= J2 fni^^ty'^' ^ (13) 

n=—OD 

with fo = 1, fn = fin- ^s ^^^^ ^^ ^'i- (H)) ^^e effcct of the noise term in Eq. (1) is 
to introduce diffusion in the phase angle 6 whose strength is characterized by the phase 
diffusion coefficient ^. 

In Sec. Ill we will solve Eqs. (11) and (12) for the case dp <^ 1, thus demonstrating the 
echo phenomenon as described in Sec. I. In Sec. IV we will present numerical solutions of 
Eq. (1) for large A^. 

III. ANALYSIS 

A. Amplitude expansion 

In order to proceed analytically we use a small amplitude expansion and obtain results to 
second order (i.e., up to quadratic in the small amplitude). This will be sufficient to obtain 
the echo phenomenon. We introduce a formal expansion parameter e, as follows, 

/ = 1 + ef^'^ + eV^') + Oie') ; (14) 



dp = edp for p = 0, 1; R = ei?(^) + e'^R^^^ + 0{e^); R^"'^ = J gfi^^du; where /("^) = 
T,nfn exp{in6). (Although we formally take e <^ 1, when we finally get our answers, the 
results will apply for e = 1 and dp = dp, ii dp <ti I.) 

B. Order e 

In linear order (i.e., 0{e)), by multiplying Eq. (11) by exp{—i6)d6 and integrating over 
9, we have for the component of f^^^ varying as e*^, 

^- + (^uJ + Ofi'^ = ^R^'^" + ^h,m , R^'^*it) = J f^'^gdco , (15) 

where /{ {uj,t) = for t < and R^^'^* is the complex conjugate of R^^\ Due to the delta 
function term on the right hand side of Eq. (15), ihidoS{t), at the instant just after the first 
delta function (denoted t = 0+), /} jumps from zero just before the delta function (denoted 
t = 0^) to the value /} {cu, 0"^) = ihido. Making use of this observation, in Appendix I we 
solve Eq. (15) for < t < r, with the result that, for K < Kc, 

/} {uj,t) = A{uj)e~^^'^~^^'^ + (a more rapidly exponentially decaying component) , (16) 

where 

A{uj) = thdo/D[-{tuj + 0] , (17) 

Dis)^l-^r-^^,forReis)>0, (18) 

and D{s) for Re{s) < is defined from Eq.(18) by analytic continuation. Since Eq. (16) 
applies for < i < r, we have that just before the application of the second delta function 
stimulus (t = r^), 

/«(a;,r-)-A(^)e"(-+«)^ (19) 

where we have neglected the second term on the right hand side of Eq. (16) on the basis 
that, due to its more rapid exponential decay, it is small compared to the first term. 

Solutions of D(s) =0 govern the stability of the state with R^^^ = 0. Let s = sq denote 
the solutions oi D(s) = with the largest real part. If -Re(so) < the state -R*-^-* = is stable, 
and a perturbation away from R^^^ = decays to zero with increasing t at the exponential 
rate Re{so). If -Re(so) > 0, then the perturbation grows and R^^^ eventually saturates into 
a sustained nonlinear state of coherent cooperative oscillatory behavior 10|, lll|. In general, 



Re{so) is an increasing function of the coupling constant K, and -Re(so) < for i^ < Kc, 
where Kc is a critical value that depends on ^ and gi^uj). Throughout this paper we shall be 
considering only the case K < Kc for which Re{so) < 0. 

It is instructive to consider ^ = 0. In that case, the first term in Eq. (16) is of con- 
stant magnitude in time, but, as time t increases, it oscillates more and more rapidly as 
a function of u. Because of this increasingly rapid variation in u, the contribution of this 
term to R^^^*{t) = J gfi duj decays in time (see Appendix I), and it does so at the same 
time-asymptotic rate as the contribution from the second more rapidly exponentially de- 
caying contribution in Eq. (16). Thus the order parameter magnitude decays away, but the 
distribution function /} can still have a component (the first term in Eq. (16)) due to the 
pulse that has not decayed away. A similar conclusion applies for ^ > provided that ^ is 
substantially less than the damping for the second term in Eq. (16). This is the source of 
the 'memory' referred to in Sec. I. It is also worth noting that the first term in Eq. (16) can 
be thought of as the manifestation of the continuous spectrum of the Kuramoto problem. 



discussed in detail in Ref . 17|] . Thus the echo phenomenon that we derive subsequently can 
be regarded as an observable macroscopic consequence of the continuous spectrum, where by 
'macroscopic' we mean that the effect can be seen through monitoring of the order parameter 
without the necessity of other more detailed knowledge of the distribution function. 
It is also of interest to consider /„ for n > 2. From Eq. (11) we obtain for \n\ > 2 

J^ + (,ncu + n^O/^) = inK/\{t) , (20) 

which does not have any contribution from the order parameter, R. For r > t > 0, Eq. (20) 
yields 

fl^'{LO,t)=inhndQe:!cg[—{inLO + n^^)t] , (21) 

for < t < r, which, similar to the first term on the right hand side of Eq. (16), also 
oscillates increasingly more rapidly with uj as t increases. At time t = t^ Eq. (21) yields 

fi^\uj, T~) = inhndo exp[-{inuj + n^O^] > (22) 

for \n\ > 2. 



C. Order e^ 

Now proceeding to 0(e'^) and again (as done in obtaining Eq. (15)) taking the e*^ com- 
ponent of Eq. (11), we have 



dt 



(_ n=—oo 



where R^^''^^*{t) = J_ °° g(uj)fl ' {uj,t)duj. The above equation is hnear in /} and is driven 
by several inhomogeneous terms appearing on the right hand side of Eq. (23) that are 

(2) 

quadratic in first order quantities. Since we are interested in the components of /^ that 
result in echoes, and since, by our previous discussion, we expect that the echoes depend 
on the presence of both stimulus delta functions (i.e., the delta function 6{t) of strength do 
and the delta function 5{t — r) of strength di), we are interested in the component of /} 

(2) 

that is proportional to the product dodi for t > t. We denote this component fl^, where 

(2) 

the subscript e stands for 'echo'. From Eq. (23) we see that for t > r, the fl^ component 

(2) 

of fl satisfied the following initial value problem 

-^ + i^u^ + OfS-^KR(^^* = 0, (24) 

+0O 



/g(.;, r+) = td, J2 h-in-i)fi'\uj, T-) , (25) 

n=—OD 

n+CO 

Rf>{t)= g{uj)f^}{uj,t)du; . (26) 



Since /« {io,T ) is proportional to rfg (see Eqs. (19) and (22)), we see that the solution of 

(2) (2) 

Eqs. (24)-(26) for fl j and -Re will indeed be proportional to dodi as desired. 
We solve Eqs. (24)-(26) by taking Laplace transforms. 



?(2)/, , ,\ _ / ^-stA2) 



fm^.s)= e-^7iJ(o;,t)rft, (27) 



RfJis)^j e-^*i?P*(t)rft , (28) 

in terms of which we obtain from Eq. (24) 

/£.(.. .,.^2'^^ + %^. (29) 

S + t, + iUJ S + t, + lUJ 



Multiplying Eq. (29) by g{uj)duj and integrating from uj = — oo to uj = +00, then yields 

i?(2)(s) = I— / ^'''^^' / g{uj)duj . (30) 



1 p+ioo+ri 

Rf>{t) = — e^'RfJ{s)ds,v>0. (31) 

^TT^ J-ioo+r? 



To find i?e (t) we take the inverse Laplace transform, 

-ioo+J7 

For the purposes of evaluating the integral (31), we recall that D{s) = has roots whose real 
parts correspond to the exponential decay rate of a response to an initial stimulus toward 
the R = state. Thus, as before in our discussion of the linear response (see Eq. (16)), 
any poles at the roots of D{s) = give contributions that we assume decay substantially 
faster with increasing t > t than the diffusion induced exponential decay rate C,- Since we 
are interested in echoes that we will find occur for t = 2r, 3r, . . ., we neglect contributions 
to Eq. (31) from such poles. Thus it suffices to consider only the contribution to Eq. (31) 
from the pole at s + C, + iu = 0. Hence Eqs. (30) and (31) yield 

RfHt) = J ^ ^"'^^''"^^:^r|^^(^)^- • (32) 

D. Echoes 

In order to see how Eq. (32) results in echoes, we recall our previous results, Eqs. (25), 
(19) and (21) for /|g(a;,r+), and combine them to obtain 

/S(^> ^'^) = dodih2hl — -— r - dodi V nKhl_^ exp[-(mcj + n^Or] , (33) 

L ^ ^'^ |n|>2 

where we have used /iq = 0, /i„ = /i*, /i^ = /} , and the first term on the right side of 
Eq. (33) corresponds to n = — 1 in Eq. (25). Putting Eq. (33) into Eq. (32), we see that we 
have an integral of a sum over terms with exponential time variations of the form 

exp{-icj[t - (1 - n)r]} exp{-^[t + {n^ - l)r} . (34) 

Considering the first exponential in Eq. (34), we see that, for large values of |t — (1 — n)r|, 
there is rapid oscillation of the integrand with a;, and the integral can therefore be expected 
to be near zero. However, such rapid oscillation is absent near the times t = (1 — n)r, at 
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which a large value of Re will occur. Since t > r, the relevant times occur for n < — 1; 
e.g., for n = —1, we get an echo at t = 2r; for n = —2, we get an echo at t = 3r; etc. 
Therefore, we henceforth replace the summation over \n\ > 2 in Eq. (33) by a summation 
from n = — oo to n = —2. 



E. Evaluation for Lorentzian frequency distribution functions 

We now consider the case of a Lorentzian frequency distribution, 

^^") = ^^^") " ^ (c-ct^+A^ = 2^ {:7^jh^) - ^:;:^jhjA)] ■ ^^^^ 

The right-most expression for gii^) makes clear that, when the previously real variable u 
is analytically continued into the complex plane, the function gii'-^) results from the sum of 
two pole contributions, one dXu = cu+iA, and one dX u = u! — iA. The quantity ui represents 
the average frequency of the distribution, while A represents the width of the distribution. 
Consideration of the Lorentzian will be particularly useful to us because the integral (32) 
can be explicitly evaluated, and also because our numerical experiments in Sec. IV will be 
for the case of a Lorentzian frequency distribution function. 

As a first illustration we consider the n = —1 term which results in an echo at t = 2r. 
We first evaluate D{s) by inserting the pole-form for gii^) into Eq. (18) and closing the 
integration path with a large semicircle of radius approaching infinity. This yields a single 
residue contribution to D{s), 

Dis) = l-^[s + ^ + z{u- zA)]-' . (36) 

Note that the solution of D{s) = occurs at 

s = -iu-U + A-^^ . (37) 

According to our previous assumptions, we require K < Kc = 2{A + C,) so that the R = 
state is stable, and (A — K/2)r — ^t ^ 1 so that we can neglect contributions from the pole 
at the root D{s) = in our approximation of (31) by (32). Using Eq. (36) and the n = —1 

(2) 

contribution to fl^ (i.e., the first term in (33)) in Eq. (32) we obtain for the echo term at 
t = 2r (denoted i?£^*(e)). 



?(2)v,^_9.-/,*/,„^„^,A / ^ expH^(t-2r)-et] 

2m [{u;-co)-tiA-fmuj-co)+tiA-f)] 



Rf>{t) = 2^hlh,d,d,A ^ . - — _, rrT.' ':.... k.. ■ m 
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For t > 2t (t < 2r) the integrand exponentially approaches zero as Im{uj) -^ — oo {Im{uj) — > 
+00), and we can therefore close the integration path with a large semicircle in the lower half 
cu-plane (upper half cj-plane). Thus the integral (38) is evaluated from the pole enclosed by 
the resulting path [i.e., the pole uj = cjq — "^(A — -j) for t > 2r, and the pole u = cJo + ^(A — y) 
fort <2r], 

R?J*{t) = M^H^^e-^(*-2^)-«*e-(^-f )l*-^-l . (39) 

From Eq. (39) we see that we obtain an echo that is approximately symmetric in shape 
about t = 2t (i.e., the envelope exp[-(A - K/2)\t - 2r|]) for ^ < (A - ^K). 

We can similarly evaluate the contribution Rmrit) of echoes at t = mr for m = 3, 4, . . .. 
For example, the result for the echo at t = 3r is 

^if = 2^2Mo^lA ^-^(3.+t)^-.^ft-3.)^(^ _ 3^) ^ (40) 

A — (ii/4) 

{exp[A(t - 3r)l , for t < 3r , 

exp -[(A - \K){t - 3r)] , for t > 3r . 
Thus, in the case ^ = 0, the shape of the pulse envelope E{t — 3r) is asymmetric about 
t = 3r, increasing at a more rapid exponential rate (namely. A) as t increases toward 3r, 
than the slower exponential rate of decrease (namely, A — {K/2)) as t increases away from 
3r. This is in contrast to the symmetrically shaped envelope exp[— (A— |i^)|t — 2r|] for the 
echo at t = 2r. 

In Appendix II we present an evaluation of R2t- (t) for the case of a Gaussian frequency 
distribution function, 

g{u)=gG{u;) ^ [2nA']-'/'exp[-{u;-uf/{2A')] . 



F. The small coupling limit 

We now consider a general frequency distribution function g{uj) but for the case where 
the coupling between oscillators is small. That is, K <^ A, where A denotes the frequency 
width of g{uj) about its mean value u = ui. In this case a good approximation is provided 
by setting K = 0. Thus D[-{iLO + 0] - 1 and Eq. (33) yields 

CO 

/i^e (u;, r+) = dodi ^ n/i*/i„+i exp[-(-mu; + n^O^] , (42) 



n=l 
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where we have replaced n by — n and used /i„ = /il„. Inserting Eq. (42) into Eq. (32) we 
obtain 

oo 

Rf^*it) = Y.^n - l)d,d,K^,K~g{t - nr)e-[("'-i)-+*]« , (43) 

n=2 

where g{t) is defined by 

/ + 00 
due-^-'g{u) . (44) 

-oo 

Thus, for K <^ A, the shape of the echoes at t = 2r, 3r, ... is directly given by the Fourier 
transform (44) of the frequency distribution function g{uj). Another point is that with 
i^ ^ 0, Eq. (1) shows that the oscillators do not interact, and the nonlinearity needed to 
produce the echo phenomenon comes entirely from the stimulus function h{6). 

IV. SIMULATIONS 

We have performed direct numerical simulations of the system (1) with a Lorentzian 
oscillator distribution (see Eq. (35)), c<) = 0, A = 1 (corresponding to Kc = 2), do = di, 
K = 1, T = 50, and ^ = 0. At t = 0" we initialize each phase 6i for i = 1,2, . . . ,N randomly 
and independently with a uniform distribution in the interval {0,2tt). We then apply the 
mapping given by Eq. (7) with dp = do to each 6i in order to simulate the effect of the delta 
function at t = 0. Next we integrate Eq. (1) for each i = 1,2,. . . ,N forward in time to 
t = T~ , again apply the mapping Eq. (7) (but now with dp = di), and we then continue the 
integration. At each time step we also calculate R{t) using Eq. (8). Figure 2 shows results 
for do = di = 1/4, and 

h{9) = sm9 + sm29 , 

for several different system sizes, A^ = 10^, 10^, and 10*^. Figure 2(a-c) shows \R{t)\ versus 
t for < t < 125. The responses to the delta functions at t = and r, as well as the echo 
at time t = 2t are clearly illustrated. The effect of lower A^ is to increase the fluctuations 
making the echo somewhat less distinct. We do not see any echo at t = 3r. This is 
in agreement with Eq. (40), since h^ = for the h{6) employed in these computations. 
Figure 3 shows a blow-up of the numerically computed echo around the time t = 2t for 
A^ = 10®, 10^, and 10'^. Also, plotted in Fig. 3 as asterisks is the result from our theoretical 
calculation Eq. (39). Reasonable agreement between the theoretical and computed echo 
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FIG. 2: \R{t)\ versus t for (a) iV = 10^, (b) iV = 10^ (c) iV = 10'', and (d) iV = 10^ showing the 
echo at t = 2t and the increase of fluctuations at lower A^. 

shapes is obtained, although the agreement is somewhat obscured by fluctuation effects 
at the smaller system sizes (A^). While our choice do = di = 1/4 might be regarded as 
questionable for applicability of the small amplitude approximation {dp <^ 1, for p = 0, 1) 
employed by Eq. (7) and by our theory of Sec. Ill, we have nonetheless evidently obtained 
good agreement between the theory and numerical experiment. Figure 4 illustrates the 
effect of varying the driving amplitude for a network of size N = 10^. For do = di = 1/8 
(Fig. 4(a)) the echo is swamped by the noise and is not seen. For do = di = 1/4 (Fig. 4(b), 
same as 2(a)) the echo seems to have appeared, but because of the noise, this conclusion is 
somewhat questionable. Finally, at the larger driving oi do = di = 1/2, the echo is clearly 
present. 

Figures 5(a) and 5(b) show the effect of changing h{9). In particular. Fig. 5(a) shows 
numerical results for do = di = 1/4, A^ = 10^, and h{6) = sin^, with all other parameters 
the same as before. Since /i2 is now zero, Eq. (39) now predicts that there is no echo, in 
agreement with Fig. 5(a). Figure 5(b) shows numerical results for do = di = 1/4, A^ = 10^, 
and 

h{e) = sin e + sin 26 + sin 3^, 

with all other parameters the same as before. Since hi, /i2 and h^ are all nonzero, Eqs. (39) 
and (40) now predict echoes at both t = 2t and at t = 3r, and this is confirmed by Fig. 5(b). 
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FIG. 3: \R{t)\ versus t blown up around t = 2t = 200, for iV = 10^ 10^, 10^, 10^ (solid curves) 
showing the increase of fluctuations at lower N. The dotted curve is the theoretical result from 
Eq. (39) with ^ = 0. 

Finally, we note that similar numerical experiments to all of the above have been repeated 
using a Gaussian g{uj), and these yield similar results (not shown). 

V. DISCUSSION 

Echo phenomena as used for MRI provide a powerful medical diagnostic tool. Echoes 
in plasmas have also been used as a basis for measuring velocity space diffusion of plasma 
particles [is!]. Thus it is of interest to consider whether there are potential diagnostic mea- 
surement uses of echoes in the context of situations that can be described by the Kuramoto 
model and its variants. For example, we note that the amplitude of the echo varies expo- 
nentially with ^, providing a possible means of determining the phase diffusion coefficient 
^. For example, the amplitude of the echo at t = 2r varies as e~^'^. Thus the log of the 
ratio of measurements of the echo amplitude using two different values of r, divided by the 
difference in the r values, provides a potential means of estimating ^. Also, as indicated by 
Eq. (43), if one can lower the coupling K sufficiently, then echoes provide a potential way of 
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FIG. 4: Simulation of 10^ oscillators for r = 100, do = di = 1/3, K = 1 = ^K^, h{9) = sin 6'. In 
this case, no echo at i = 2r = 200 is observed. 



determining the oscillator frequency distribution function g{(jj). In particular, for low K the 
distribution g{uj) is directly given by the inverse Fourier transform of the echo profile. On 
the other hand, we have seen from the simulations in Sec. IV that finite A^ leads to noise-like 
behavior that may compromise such attempts. We also note that the Kuramoto model is an 
idealization, and application to any given situation may require modifications of the model 
and theory to more closely correspond to the situation at hand. We, nevertheless, feel that 
consideration of echoes for diagnostics may be of potential use. 

Furthermore, these phenomena are of theoretical interest from at least two points of 
view. First, as mentioned in Sec. Illb, the memory required by the echo phenomenon can 
be thought of as leading to a macroscopically observable consequence of the continuous 
spectrum 171] of the Kuramoto model. A second point of theoretical interest relates to the 



recent work in Ref. 



19( 1 . In that paper it was shown for a general class of initial conditions 



that are on a certain manifold of the infinite dimensional state space of the Kuramoto system, 
that the future time evolution of the order parameter is determined by the current value 
of the order parameter. In particular, there is an ordinary differential equation describing 
the order parameter evolution. The echo phenomenon provides an example showing that, if 
initial conditions do not lie on the specified manifold of Ref. [l9|, other behavior can occur. 
In particular, well after the second stimulus (at t = r) and well before the occurrence of the 
first echo (at t = 2r), the order parameter is essentially zero, yet it does not remain zero 
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FIG. 5: Simulation of 10^ oscillators for r = 100, do = di =??, K = I = \Kc, h{e) = sin6i + 
sin 2^ + sin 3^. In this case, echoes are seen at t = 2t = 200 and at t = 3r = 300. The inset shows 
a blow-up of the numerical result for the echo shape at t = 3r with the theoretical result, Eq. (41), 
superposed (dotted curve). 



as would be predicted for initial conditions on the manifold of Ref. [19| for K < Kc- This is 
discussed further in Appendix III. 

In conclusion, we hope that our work will stimulate experimental groups to investigate 
the type of situations we have addressed. 

This work was supported by ONR (N00014-07-1-0734) and NSF (PHY0456249). 



Appendix I: Linear Analysis 

In this Appendix we solve Eq. (15) for < t < r to obtain the solution (16) and (17) for 
K < Kc- Taking the Laplace transform, u{s) = J^ u(t)e~'^^dt, Eq. (15) yields 



?(i) 



/r(^,^) 



^R^^\s)+zhdo)/{s + ^ + zu;) 



(45) 



where Rl (s) denotes the Laplace transform of R^^^*{t). Multiplying Eq. (45) by g{uj)di 



.LO 



and integrating from uj = —oo to uj = +oo, we obtain 



R^}\s) = ih^doI{s)/D{s) 



(46) 
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where I{s) = j_^ duog{uo)/{s + i + ioo) and D{s) = 1 — {K/2)I{s). Inserting Eq. (46) in 
(45) gives 

f['\uj,s) = ihrd,[D{s){s + i + ioo)]'^ . (47) 

As noted in Sec. Illb, /} (a;, s) has poles in s at the zeros of D{s) and at s = —{iuj + ^). 
These yield time dependences of the inverse Laplace transform of /} (see Eq. (27)) that 
vary as e'^^^ and as e"'-*'^^^-'*, respectively, where Sq denotes the root of D{s) = with the 
least negative real part. For t ^ r and — [i?e(so) + ^]t ^ Ij we can neglect the contributions 
from poles arising from roots of D{s) = 0, and use only the contribution from the pole at 
s = —{iuj + ^). From Eqs. (27) and (47) this yields 

if\u,t)^ih,d,e-^'-+^^'/D[-{iu + i)] , (48) 

thus confirming Eqs. (16) and (17). 

Appendix II: Echo at t = 2r for Gaussian g{uj) 

We consider the case g{uj) = gci^) = (27rA^)"^/^ exp[-(cj - cj)^/(2A^)]. Putting this 
expression for g{uj) and the n = —1 contribution to flJ. (i.e., the first term in Eq. (33)) into 
Eq. (32) we have, 

^,W. h,hM r^ exp-{ i^-H-t:^-)1- +f(t-2r)^ + ..(t-2.)-et} 

(49) 
The collective damping rate is determined by the root of D{s) = with the least negative 
real part. Denote this root s = sq where 

So = -{iuJ + ^ + 7o) = -(^^0 + , ^0 = ^ - ^7o , (50) 

where 70 > is real. Letting F{uj) = D[—{iuj + ^)], continuing this function from real uj into 
the complex cu-plane, and expanding around u; = o^o, we have 

F{lu) = {lu- iUo)v + 0[iLU - cuo)'] , (51) 

where 77 is a complex constant. Letting F^:{uj) denote the continuation of the function of 
the real variable u in Eq. (49), D*[—{iuj + ^)], into the complex u;-plane, we have that this 
function has a zero at u; = cUq, 

F^iu) = {uj- u;)r]* + 0[{uj - u;*f] . (52) 
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Considering the oscillatory uj variation in the numerator of the integrand of Eq. (49) to be 
rapid (valid for K <^ 70), we can approximate the integral by the saddle point method, 
where the saddle point is at 

ujsp = oj — iA'^{t — 2t) , 

and the steepest descent path through u = uogp runs along the horizontal line Im{uj) = 
— A^(i — 2r) from Re{uj) = —00 to Re{uj) = +cxd (see Fig. 6). From Fig. 6(a) we see that for 
A^|t — 2r| < 7o, the poles at u = lj ± i'jQ are not intercepted by the steepest descent path, 
while for A^|t — 2r| > 70 one of the poles is intercepted (e.g.. Fig. 6(b)). In the case where 
a pole is intercepted, its contribution dominates the contribution from the saddle point by 
virtue of its time dependence, e"''''''*"^'^', as opposed to the saddle point contribution time 
dependence, e~2^ (*-2t) _ Thus we obtain 

(2)*. . -u 9 w. f e-^(*-2-)' for |t-2r| < 270/A2 , ^ ^ 

I e~^ol*-2^l for |t - 2r| > 270/A2 . 

Near 70 = A^|t— 2r|/2, the pole is near the saddle point, and a uniform asymptotic expansion 
of the integral (49) is necessary to obtain the transition between the two forms in Eq. (53). 



Appendix III: Further Discussion of Ref. 19|] 



In Ref. [l9l], a broad class of noiseless (e.g., ^ = in Eqs. (3) and (11)) globally coupled 
systems of phase oscillators was studied. The simplest example of this class is the Kuramoto 
model. Reference [l9| considered Lorentzian g{uj) and a special class of initial conditions. 
Referring to Eq. (12), these initial conditions are of the form, 

/„(a;,0) = a"(u;) , for n > , (54) 

and /rt(^, 0) = /*„(ci;, 0) , for n < , where |a(u;) | < 1 for cu on the real axis, a{uj) is analytic 
in \m.{uo) < 0, and |tt(u;)| — i> as Im(c<j) — > —00. Under these conditions, Ref. 19| shows that 
the order parameters (or parameter), see Eq. (12), that describe the nonlinear, macroscopic 
time evolution of the given system satisfy a finite set of ordinary differential equations 
in time. Thus the order parameter dynamics is low dimensional, while the dynamics of 
the full system determining the evolution of the distribution function f{uj,6,t) is infinite 
dimensional 19|. For example, for the Kuramoto problem with the above conditions satisfied, 
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FIG. 6: (a) Steepest descent path (dashed) through the saddle point uj = ujgp for A^|t — 2r| < 7*. 
(b) The steepest descent path (dashed) for A^(t — 2r) < —7*. The dominant poles at the roots 
of F{uj) = and F*(u;) = are shown as crosses, where in (b) the steepest descent path has 
intercepted the pole oo = ujq resulting in a pole contribution to figr (0- 



Ref. \m shows that 



dR/dt + ( A - ]-k\ R + iif A|i?pi? 







(55) 



where we have taken a; = in Eq. (35). 

A consequence of Eq. (55) is that for K < 2A = Kc, \R{t)\ decreases monotonically to 

a.UQwed in the eeho-^i^ienQmena we-discuss-in^he-p-resent-pa 
In particular, in Fig. 1, |-R(t)| is small between t = t and t = 2r, but then increased to form 
the echo in the vicinity of time t = 2r. Referring to Eq. (34) and our subsequent discussion, 
we see that this is because there is a component of fi{uj,t) that varies as exp[—iuj{t — 2t)]. 
Identifying fi{uj,0) in the linear problem with a{uj) in the nonlinear problem [Eq. (54)] 
and considering to as a new initial time (shift time so that to goes to t = 0), we see that 
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FIG. 7: (a) Steepest descent path (dashed) through the saddle point co = ujgp for A^|t — 2r| < 70. 
(b) The steepest descent path (dashed) for A^(f — 2r) < —70- The dominant poles at the roots 
of F{u)) = and -F*(u;) = are shown as crosses, where in (b) the steepest descent path has 
intercepted the pole a; = ojo resulting in a pole contribution to -Rgr (^)- 

a{uj) ~ exp[— ^^^(to — 2r)]. If we take to to be such that t < to < 2t and |-R(to)| is small, then 
a{u!) does not satisfy the condition of Ref. 19| that a{uj) ^ as 1171(00) -^ —00. However, 



if to > 2t, then it does. Thus the increase of |-R(t)| occurs only when the hypothesis under 
which Eq. (55) was derived does not hold. 

More generally, consider an initial condition for the original Kuranioto problem (without 
stimuli or noise) where fi{uj, 0) isanalytic on the real a;-axis. Expressing fi{uj, 0) as a Fourier 
integral transform, we have 



/i(^,0) 



iLorj 



k{r])dr] 



(56) 



where k{ri) is the Fourier transform of fi{uj, 0). Since fi{uj, 0) is analytic in u, k{ri) decreases 
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exponentially for sufficiently large rj, 

\k{rj)\ < He-f'^ if r]>rjo, (57) 

for some set of positive constants H, f3, rjQ. Using the Laplace transform technique (as in Ap- 
pendix I), it can be shown that the solution to the linearized initial value Kuramoto problem 
contains a component of fi{uj, t) of the form exp{—iujt)fi{uj, 0), which we can express using 
Eq. (56) as 

/t /"OO 

e-'-^(t-^)k{r^)dr^+ e''^^'^-'^k{r])d7] . (58) 

OO Jt 

Setting t = to and regarding t = to as a new initial condition time, we note that the initial 
condition consists of two terms, namely the ffist and second integrals on the right hand side 
of Eq. (58). For to > ^70 sufficiently large, the second integral is smaller than the ffist by a 
factor of order exp(— /5to). Furthermore, the ffist integral satisfies the condition fi{uj, to) -^ 
as Im(u;) —>■ — oo [because {t] — to) > for the ffist integral], while the second integral does 
not. Thus, if we choose to shift what we designate as the initial time to sufficiently large to, 
then aside from an exponentially small component of order exp(— /3to), the initial condition 
obeys the requirement of Ref. [l9|] that fi{uj,to) goes to zero as Im(u;) — > — oo. 
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